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O Abstract 

The Mars Atmosphere-Ice Coupler MAIC-2 is a simple, latitudinal model, which 
consists of a set of parameterisations for the surface temperature, the atmospheric 
water transport and the surface mass balance (condensation minus evaporation) of 
water ice. It is driven directly by the orbital parameters obliquity, eccentricity and 
solar longitude (L s ) of perihelion. Surface temperature is described by the Local 
Insolation Temperature (LIT) scheme, which uses a daily and latitude-dependent 
radiation balance. The evaporation rate of water is calculated by an expression for 
CN free convection, driven by density differences between water vapor and ambient air, 

^T) the condensation rate follows from the assumption that any water vapour which ex- 

ceeds the local saturation pressure condenses instantly, and atmospheric transport of 
water vapour is approximated by instantaneous mixing. Glacial flow of ice deposits 
is neglected. Simulations with constant orbital parameters show that low obliquities 
• i-h favour deposition of ice in high latitudes and vice versa. A transient scenario driven 

by a computed history of orbital parameters over the last 10 million years produces 
essentially monotonically growing polar ice deposits during the most recent 4 million 
years, and a very good agreement with the observed present-day polar layered de- 
posits. The thick polar deposits sometimes continue in thin ice deposits which extend 
far into the mid latitudes, which confirms the idea of "ice ages" at high obliquity. 
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1 Introduction 



On time scales of 10 5 -10 6 years, Mars has experienced large periodic changes of the orbital 
elements obliquity, eccentricity and equinox precession. These changes have an impact on 
the Martian climate. The obliquity determines the strength of the seasons and the latitu- 
dinal distribution of mean solar insolation. The eccentricity determines the magnitude of 
the asymmetry of insolation with season, and the equinox precession determines the timing 
of the asymmetry of solar insolation with season. On Earth, the so-called Milankovitch 
cycles of much weaker orbital changes with periods of 20, 40 and 100 ka are considered 
driving forces for climate variations like the glacial/interglacial cycles. It can, therefore, be 
expected that the main Martian ±10° obliquity cycles with periods of 125 ka and 1.3 Ma 
and the secular shift from high to low average obliquities at 4-5 Ma ago ( Laskar et al.|2004 



shown in Fig. [I]) have significant impacts on the climate and the polar layered deposits 
(PLDs) due to large insolation changes in the polar regions. 
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Figure 1: Martian obliquity for the last 10 Ma (Laskar et al. 2004). 



This idea is supported by the presence of light-dark layers in the PLDs, which are 
exposed in the surface troughs and close to the margins, and which are actually the reason 
for the term "polar layered deposits" . These layers indicate a strongly varying dust content 
of the ice due to varying climatic conditions in the past. During periods of high obliquities, 
insolation in the polar regions is large, which entails higher sublimation rates of superficial 
ice of the PLDs and probably of permafrost in the ground. This may lead to the formation 
of a thicker and dustier atmosphere, so that dust accumulates on the PLDs. By contrast, 
during periods of low obliquities, the atmosphere is thin and dust deposition is low, so 



that clean ice forms at the surface of the PLDs. Along this line of reasoning, Head et al. 



(2003) presented evidence for past glaciation in the mid-latitudes and suggested that Mars 
experienced "ice ages" during periods of high obliquity like that from about 2.1 to 0.4 Ma 
ago (with obliquity maxima of ~ 35°). These ice ages were supposedly characterised 
by warmer polar climates, enhanced mass loss of the PLDs due to sublimation and the 
formation of metres-thick ice deposits equatorward to approximately 30°N/S. 

In a number of studies, General Circulation Models (GCMs) have been applied to 
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). These models. 



all derivatives of Earth GCMs, solve the equations of fluid dynamics and thermodynamics 
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and include e.g. the processes of radiative transfer, cloud formation, regolith-atmosphere 
water exchange, and advective transport of dust and trace gases. However, they have 
essentially been designed to simulate the present-day atmosphere in as much detail as pos- 
sible, and thus are computationally too expensive to permit long-term paleoclimate studies. 



Segschneider et al. (2005) and Stenzel et al. (2007) adapted an Earth System Model of 



Intermediate Complexity (EMIC) to Mars. This "Planet Simulator Mars" (PlaSim-Mars, 
formerly called "Mars Climate Simulator" ) allows in principle simulations over longer, cli- 
matological time scales. So far, only scenarios for present-day conditions and varied obliq- 
uity angles have been considered, and the impact on the PLDs has been studied by cou- 
pling PlaSim-Mars with the three-dimensional, dynamic/thermodynamic ice-sheet model 
SICOPOLIS (http://sicopolis.greveweb .net/[ ). In addition to that, simple one-dimensional 
models have been used to study specific processes that do not require full solution of the 
dynamic equations. Examples are the radiative transfer model by |Gierasch and Goody 



(1968), the energy balance model by Armstrong et al. (2004), regolith-atmosphere water 
exchange (Jakosky 1985), and formation of water ice clouds (Michelangeli et al. 1993). 

In this study we aim at simulating the surface glaciation of the entire planet with 
a simple model that depends only on latitude and time. This model, termed the Mars 
Atmosphere-Ice Coupler Version 2, or MAIC-2 in short, consists of a set of parameterisa- 
tions for the surface temperature, the atmospheric water transport and the surface mass 
balance (condensation minus evaporation) of water ice. It is driven directly by the orbital 
parameters obliquity, eccentricity and solar longitude (L s ) of perihelion, which were pub- 
lished by Laskar et al. (2004) for the period from 20 million years ago until 10 million years 
into the future. MAIC-2 is applied to two different kinds of scenarios, namely (i) scenarios 
with orbital parameters kept constant over time, and (ii) transient scenarios forced by the 
history of orbital parameters over the last 10 million years. The evolution of surface glacia- 
tion is studied for these scenarios under the simplifying assumption of negligible glacial 
flow, so that changes of local ice thickness are exclusively determined by the local surface 
mass balance provided by MAIC-2. 



2 Model MAIC-2 

The design of MAIC-2 is illustrated schematically in Fig. [2] All quantities are latitude 
(<£>) and time (t) dependent, with the exception of the atmospheric water content. Since 
instantaneous mixing is assumed, only the (time dependent) global mean water content is 
modelled. The formulation of the different processes is detailed in the following sections. 

2.1 Timekeeping and orbital position 

We use the timekeeping of PlaSim-Mars (see Sect. [I]), in which a Martian year consists 
of 12 months of 56 days (sols) length. A day has 24 hours of 61.5 minutes length, and a 
minute is 60 SI seconds long. 

In order to compute the relation between solar longitude L s and time t, let us identify 
the beginning of a Martian year with the northern hemisphere vernal equinox (L s = 0°). 



3 



Insolation 



C0 2 cover 



Surface temperature 



Evaporation 


Condensation 








Net mass balance 




Ice thickness 


Atmospheric water content 



Figure 2: Schematic illustration of the quantities and processes modelled by MAIC-2. 
The orbit of Mars around the sun is described by the ellipse 



1 + e cos ip 

where r is the distance Sun - Mars, ip the orbital position measured from the perihelion 
("true anomaly"), p the semi-latus rectum and e the eccentricity. By combining this 
equation and the conservation of angular momentum I, 

(where m denotes the mass of Mars) , we find 

ip = = -^(1 + ecosV) 2 = fi(l + ecos^) 2 (fi : =J_). (3) 
mr z mp* v mp AJ 

Let L s p be the solar longitude of perihelion, then 

if) = L s - L SiP , (4) 

and since L S)P varies only slowly over time, we can rewrite Eq. ([3| as 

L s = Q[l + ecos(L s - L SiP )] 2 . (5) 

This equation can be integrated numerically over a full Martian year (from vernal equinox 
to vernal equinox) for any values of e and L s p by a simple Euler-forward scheme. The 
initial condition is L s = 0°, and the parameter Q is adjusted iteratively such that after 
one Martian year the orbit is closed (L s = 360°), starting from the initial guess fi^t = 
27r/(l Martian year) [which is the correct value for a circular orbit with e — 0] . 

For present-day conditions (e = 0.093, L SjP = 251.0°), the result is shown in Fig. [3] 
Since the eccentricity is much larger for Mars than for Earth, the relation between L s 
and t is significantly different from a linear one. The deviation becomes as large as 21° 
(L s = 158.97° instead of 180° for day of year 336, leading to a lag of the northern autumnal 
equinox by 37.7 Martian days). 
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Figure 3: Relation between solar longitude L s and time t for present-day conditions (solid 
line). For comparison, the dashed line shows the linear relation for a circular orbit. 



2.2 Surface temperature 

In order to derive a parameterisation for the daily mean local surface temperature T(ip, t) 
(depending on latitude ip and time t), we start with the radiation balance 



aT 4 



A)F. 



(6) 



where a is the Stefan-Boltzmann constant (er = 5.67 x 10 8 Wm 2 K 4 ), A is the surface 
albedo (globally A = 0.3 assumed) and F is the local daily mean insolation as a function 



of the orbital parameters obliquity, eccentricity and solar longitude of perihelion (Laskar 



et al. 2004). In the absence of seasonal CO2 frost, Eq. (j6| provides reasonable results for 
the surface temperature. However, the equation does not account for the fact that at 



T = T, 



cond 



a - In P [hPa] 



(7) 



(where P is the surface pressure, a = 23.3494 and b = 3182.48 K; James et al. 1992) 



condensation of the atmospheric CO2 (formation of the seasonal ice cap) sets in. The 
seasonal variation of the surface pressure is neglected, and we use the global annual mean 
value P = 700 Pa instead. For this value, Eq. (J7l) yields a condensation temperature of 
Tcond = 148.7 K. Since the atmosphere never freezes out completely, this value constitutes 
the minimum of surface temperatures which can be realised. 

In order to find out when the seasonal CO2 ice cap at a given latitude ip is present, 
and therefore T = T con d holds, the seasonal cap is assumed to exist between the onset of 
the polar night (tdusk) and an unknown time t after the end of the polar night (tdawn)- 
During the polar night, condensation takes place, and the amount of formed CO2 frost 
corresponds to the energy (per area unit) 



W, 



cond 



*dusk 



°" T c 4 ond dt - 



After dawn, the solar insolation causes the CO2 frost to evaporate, which requires the 
energy 



W, 



evap 



t 

^dawn 



l-A)F-aT l 



cond 



dt. 



(9) 
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The time t at which the CO2 frost has evaporated completely follows from 



W, 



evap 



w, 



cond 



(10) 



The scheme defined by the radiation balance (J6|, modified by CO2 condensation follow- 
ing Eqs. (|7|)-([l~0|), is referred to as Local Insolation Temperature scheme (LIT); it was first 
laid down by B. Grieger (2004; talk at 2nd MATSUP workshop, Darmstadt, Germany). 
The resulting daily mean surface temperatures over one Martian year for present-day con- 
ditions are shown in Fig. |4| They agree well with the data given in the Mars Climate 
Database (Lewis et al. 1999[ ). The most notable discrepancy is that the LIT scheme over- 



predicts the summer temperatures at and very close to the poles. 
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Figure 4: (a) Daily mean surface temperature (in K) of the LIT scheme for present-day 



conditions, (b) Same, but from the Mars Climate Database (Lewis et al. 1999) 



The Mars Atmosphere-Ice Coupler with the LIT scheme and the simple treatment of 



the surface mass balance described by Greve et al. (2004) and Greve and Mahajan (2005) 



is referred to as "MAIC-1.5". It was used by these authors to drive simulations of the 
north polar layered deposits with the ice-sheet model SICOPOLIS. 



2.3 Saturation pressure of water vapour 

The water-vapour saturation pressure P sat is obtained from the Clausius-Clapeyron rela- 
tion, which can be integrated only approximately. Different approximations are available; 
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we use the Magnus- Teten formula for water vapour over ice (Murray 



1967) 



P sat (T) = Aexp(^^) , (11) 
with A = 610.66 Pa, B = 21.875, C = 7.65 K and T = 273. 16 K, which has also been 



implemented in the Planet Simulator Mars (Stenzel et al. 2007). 
2.4 Evaporation 



Ingersoll (1970) discussed the water vapour mass flux in the Martian carbon dioxide at- 
mosphere. The evaporation rate E of water from the surface, expressed as a mass flux in 
kg m~ 2 s _1 , is 

£ = £ x0.17A IfP fl(^M^y /3 , (12) 



V 



where Eq is the evaporation factor (default value equal to unity), Arj the concentration 
difference at the surface and at distance, p the atmospheric density, Ap the CO2 density 
difference at the surface and at distance, D the diffusion coefficient of water in CO2, g the 



acceleration due to gravity and v the kinematic viscosity of C0 2 (cf. also Sears and Moore 



2005). The term A77 is given by 



„sat T\,f p 
p^_ = 

1 p M C P ' v ; 

where p 1 ^ 1 is the saturation density of water vapour at the surface temperature T and M w 
and M c are the molecular weights of water and carbon dioxide, respectively. The terms p 
and Ap/p are calculated by applying the ideal gas law, 

= MeP Ap = (M c - M w ) P sat 
9 RT ' p M c P-{M c -M w )P sat ' 1 j 

where R is the universal gas constant. Parameter values are given in Table [TJ 

Quantity Value 

Gravity acceleration, g 3.72 ms -2 

Diffusion coefficient, D 1.4 x 10 _3 m 2 s _1 

Kinematic viscosity of CO2, v 6.93 x 10~ 4 m 2 s _1 

Universal gas constant, R 8.314 J mol -1 K _1 

Molar mass of water, M w 1.802 x lO^kgmor 1 

Molar mass of C0 2 , M c 4.401 x 10" 2 kgmor 1 



Table 1: Physical parameters for the evaporation model of MAIC-2. 



Sears and Moore (2005) state that the evaporation rate of ice is probably about half 
that of liquid water. In addition, any significant evaporation of the dusty ice of the PLDs 
will lead to an enrichment of dust at the surface, thus producing an insulating layer which 
hampers further evaporation. These effects can be accounted for by setting the evaporation 
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factor E in Eq. ( 12 ) to a value less than unity, and we will use a standard value of E = 0.1 



(this choice will be detailed below in Sect. 3.2) 



The dependence of evaporation on surface temperature and pressure is illustrated in 
Fig. [5j The temperature dependence is evidently very strong, while the pressure depen- 
dence is weak. Owing to the strong temperature dependence and the short reaction time 
of evaporation on changing conditions, it is not sufficient to calculate evaporation rates on 
the basis of daily mean temperatures. Therefore, we parameterise the daily cycle Tbc as 
follows, 

T DC (^, t) = T(<p, t) + f(<p) cos (^-) , (15) 



lsol 



with the amplitude 



T(<p) = r, 



EQ 



V90° 



(16) 



The amplitude at the equator is set to Teq = 30 K. This choice provides a good agreement 
with the amplitudes measured by the surface missions Mars Pathfinder (19° N, T ~ 30 K) , 
Viking Lander 1 (22°N, f ~ 30 K) and Viking Lander 2 (48°N, f ~ 25 K) ( |Tillman||200l| ) 
as well as the requirement T = K for the poles (Fig. ^ . 

For high temperatures (e.g., Tdc > 272.8 K for P = 700 Pa), due to the increasing 



saturation pressure P sat , the term Ap/ p of Eq. ( 14 ) becomes larger than unity, goes through 



a singularity and then becomes negative. This means that Ap/p loses its physical meaning. 
In that case, we correct the problem by resetting Ap/p to unity. 



The above equation (12) is valid for an ice cap in contact with the atmosphere. By 



contrast, for the case of ground ice, we assume that the evaporation rate is reduced with 
increasing thickness H reg of the ice-free regolith layer (which separates the atmosphere 
from the ground ice), 



E — > E x exp 



rcg 



7: 



(17) 



reg 



where 7 reg is the regolith-insulation coefficient, chosen as 7 reg = 0.1m. 
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Figure 6: Amplitude T of the daily cycle of the surface temperature according to Eq. (16) 



Comparison with the data for Mars Pathfinder (PF), Viking Lander 1 (VL1) and Viking 



Lander 2 (VL2) (Tillman 2001) as well as the north pole (NP). 



2.5 Condensation 

The water content u in the atmosphere is expressed as an area mass density in units of 
kgm -2 . Multiplied with the gravity acceleration g, this becomes equivalent to the partial 
pressure of water vapour at the surface. Thus we compare this pressure to the water vapor 
saturation pressure P sat and assume that all excess water condenses instantly, 

If gu > P sat (T) : excess water [u - P sat (T)/g] 

determines condensation rate C , (ig) 

else : (7 = 0. 



Note that this is a very simplistic approach because in reality condensation takes place 
higher in the atmosphere where the temperature may differ considerably from the surface 
temperature. 



2.6 Transport 

As mentioned above, the water content u in the atmosphere is an area mass density. Since 



the evaporation E (cf. Sect. 2.4) and condensation C (cf. Sect. 2.5) are expressed as mass 
fluxes in units of kgm -2 s -1 , its evolution is governed by 

^ = -V-F + E-C, (19) 

where F is the horizontal water flux in units of kgmT 1 s _1 . 

We approximate the atmospheric water transport by instantaneous mixing (on a time 
scale of several sols). This can formally be obtained by assuming a diffusive flux, 

F = -KVu, (20) 

with the limit of infinite diffusivity, K — >■ oo. 

The MAIC version with the LIT scheme of Sect. 12.21 and the surface mass balance of 
water ice that results from Sects. I2.3H2.6I is referred to as "MAIC-2" . 
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2.7 Ice evolution 



With the condensation C and the evaporation E, the net mass balance a net of the ice caps, 
expressed in units of ms" 1 ice equivalent, is 



C-E 



b net 



(21) 



where p ice = 910 kg m 3 is the density of ice. For a static model (glacial flow neglected), 
the evolution of the ice thickness, H, is then simply 



dH 



<^net • 



(22) 



Note that we allow for negative ice thicknesses (H < 0). Such a situation is interpreted as 
ground ice under an ice-free regolith layer of thickness H ieg = \H\. The thickness of the 
ground ice layer itself is undefined. 

The validity of the assumption of negligible glacial flow is debatable. On the one hand, 
modelled ice flow speeds on Mars are generally slow, even during periods of high obliquity 
( Greve et al.|2004 Greve and Mahajan|2005 ). One the other hand, locally enhanced glacial 
flow may occur near chasmata and troughs of the PLDs ( Hvidberg||2003 Greve|[2008| , and 
Winebrenner et al. (2008) argue that the overall topography of Gemina Lingula (also 
known as Titania Lobe), the lobe of the northern PLDs to the south of Chasma Boreale, 
was likely shaped by past glacial flow. In this study, we will stick to the simple, static ice 
model, but consider the inclusion of glacial flow for future work. 



3 Simulations 

We will now discuss the application of MAIC-2 to two different sets of scenarios. The first 
set is of "academic" nature with orbital parameters kept constant over time, whereas the 
second one employs a realistic, time-dependent forcing over the last 10 million years. In 
order to carry out these simulations, a finite-difference/finite- volume discretisation of the 
model equations of MAIC-2 has been derived (see Appendix [A] for details) and coded in 



the Fortran program maic2.F90 (available as free software at http://maic2.greveweb.net). 



Instantaneous mixing of water vapour in the atmosphere is assumed (diffusivity K — > oo), 
a time step of At = 0.02 a (« 7 sols) and an equidistant grid spacing of Aip = 1° are 
chosen (the formulation in Appendix |A| allows also for non-equidistant grid spacings), and 
the initial ice distribution is a layer of 19 m thickness on the entire surface of Mars. The 
latter value accounts for the ice inventory of the present-day PLDs, ~ 1.14 x 10 6 km 3 in 



the north (Grima et al. 2009) and ~ 1.6 x 10 6 km 3 in the south (Plaut et al. 2007). 



3.1 Constant orbital parameters 

Simulations #1-4 have been run over 10 million years with constant orbital parameters as 
follows: 

• Obliquities: = 15° (#1), 25.2° (present-day value, #2), 35° (#3) and 45° (#4). 
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• Eccentricity: e = 0.093 (present-day value, #1-4). 

• Solar longitude of perihelion: L s p = 251.0° (present-day value, #1-4). 

The evaporation factor in Eq. (12) is set to the standard value E = 0.1 for all four 
simulations. 

The resulting net mass balance of water ice in the first Martian year of simulation #2 
(all parameters at their present-day values) is shown in Fig. [7j The distribution resembles 
that of the surface temperature (Fig. f4k) . The seasonal CO2 caps are efficient cold traps for 
atmospheric water vapour, which leads to strongly positive mass balances in the Martian 
polar regions for most of the year. By contrast, in the lower latitudes negative mass 
balances prevail, so that the initial ice layer of constant thickness is redistributed towards 
the poles. 
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Figure 7: Simulation #2: Net surface mass balance (in mmiceequiv. a -1 ) in the first 
Martian year. The thick contour shows the equilibrium line (zero mass balance). 

The evolution of the ice thickness over the entire simulation time of 10 7 years for all 
simulations is depicted in Fig. [8] The obliquity of simulation #1 (0 = 15°) is approximately 
equal to the minimum value which occurred during the last 4 Ma (Fig. [I]). The resulting 
evolution of the ice thickness is shown in Fig. [8^,. The simulation produces bipolar ice 
deposits, somewhat more pronounced in the north than in the south. After 10 7 years, the 
simulated polar deposits reach maximum thicknesses of ~ 500 m, which is about a factor 
5 thinner than the observed PLDs at present. 

For simulation #2 (Fig. ^jp), it is striking that the redistribution of ice towards the poles 
is strongly antisymmetric. Until 10 4 years, MAIC-2 produces more pronounced ice deposits 
in the northern hemisphere and less pronounced ones in the southern hemisphere. At 10 5 
years and later, the ice migrates entirely to the northern hemisphere and concentrates 
around the north pole. In fact, the simulated north polar deposits at 10 7 years resemble the 
currently existing PLDs in extent and thickness (see below, Fig. [ToJd), while the simulated 
south polar region is ice-free. 

The reason for this behaviour is the hemispheric asymmetry of the daily mean surface 
temperature (Fig. [4]). As a consequence of the closer proximity of Mars to the Sun during 
southern summer, the southern summer is distinctly warmer than the northern summer. 
This leads to large evaporation rates during southern summer and thus a large amount 
of water vapour in the atmosphere, which is trapped preferably in the winter-cold high 
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Figure 8: (a) Simulation #1 (0 = 15°), (b) #2 (0 = 25.2°), (c) #3 (0 = 35°) and (d) #4 
(<b = 45°): Evolution of the ice thickness H. Note the different scales of the if-axes. 



northern latitudes. Conversely, the northern summer is less warm, evaporation rates are 
lower, and therefore the potential for water ice accumulation in the south polar region is 
much smaller (see also Fig. [7]). This holds also for simulation #1; however, the effect is 
less pronounced as a result of the weaker seasonal cycle due to the lower obliquity. Hence 
the hemispheric asymmetry is much weaker for simulation #1 than for simulation #2. 

For simulation #3, the obliquity (0 = 35°) is approximately equal to the maximum 
value during the last 4 Ma. Figure [SJ: displays the resulting evolution of the ice thickness. 
The result is similar to that of simulation #2; however, the concentration of ice around 
the north pole is more extreme. The simulated north polar deposits at 10 7 years are as 
thick as 6.7 km, almost 2.5 times thicker than their observed present-day counterparts. 

The obliquity of simulation #4 (0 = 45°) was reached in several maxima during the 
period of high average obliquity prior to 5 Ma ago. This makes the seasons very extreme. 
Like in simulations #1-3, the ice is preferentially deposited in the northern hemisphere 
(Fig. [8]i) due to the warmer southern summers. However, now the poles receive substan- 
tial insolation during the respective summer season, so that the northern hemispheric ice 
deposits are not thickest at the north pole any more. Instead, ice deposition is favoured in 
the low latitudes, and beyond 10 4 years simulation time the deposits develop a thickness 
maximum as far south as ~ 10°N. 
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3.2 Evolution over the last ten million years 



Simulations #5-8 attempt at providing realistic, time-dependent scenarios over the last 
millions of years. To this end, they have been run from 10 Ma ago until today, driven by 
the history of orbital parameters (obliquity, eccentricity, solar longitude of perihelion) by 
Laskar et al. (2004). The evaporation factor in Eq. (12) is set to the values E = 0.05 



(#5), 0.1 (standard value, #6), 0.2 (#7) and 0.3 (#8) 
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12.05 


Obs. 




1.14 x 10 6 


1.6 x 10 6 


2773 


228 


5 





Table 2: Simulations #5-8: Evaporation factor Eq, present-day volume of the north and 
south polar layered deposits (Vnpld, Vspld), present-day ice thickness at the north and 
south pole (i^NP, Hsp). The last row shows the observed volumes ( Plaut et al.|2007 , Grima 



et~aL]|2009t and ice thicknesses QZuber et al. | [l998| |Smith et al.||l999[ |Greve et al.||2004 
also the caption of Fig. 10). The misfits J have been computed according to Eq. (23). 



sec 



For the present (t — 0), all simulations produce bipolar ice deposits. An overview of the 
results is given in Table [2j It shows that, with increasing evaporation factor Eq, the volume 
of the northern deposits decreases, while the volume of the southern deposits increases. 
This holds essentially also for the ice thicknesses at the poles. The only exception is the 
north polar thickness of simulation #8 compared to #7, which decreases by ~ 6% even 
though the ice volume increases by ~ 4%. These distinctive trends allow to identify the 
most suitable value of Eq. To this end, we define the misfit J as follows, 



J 




(23) 



The standard deviations are introduced to make the various contributions to J commen- 
surate. They are computed from the four respective values of simulations #5-8, 



'-'VsplD 

°"Hnp 
a H SP 



0.134 x 10 6 km 3 
0.126 x 10 6 km 3 
301m, 
692 m. 



(24) 



The resulting misfits J are listed in the last column of Table [2] Simulation #6 with 
Eq = 0.1 produces clearly the best agreement (and thus Eq = 0.1 is used as standard value 
throughout this study). The errors of the volumes are as small as 1.6% for the NPLD and 
3.2% for the SPLD, the ice thickness is 13.3% too small at the north pole and 19.6% too 
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Figure 9: Simulation #6: (a) Maximum ice thickness (the circle marks correspond to the 



time slices shown in Fig. 10). (b) Latitude of maximum ice thickness. 



large at the south pole. Given the simplicity of the MAIC-2 model, this is a remarkably 
good overall agreement. 

In the following, the best-fit simulation #6 will be discussed. The maximum ice thick- 
ness and its position on the planet are shown in Fig. [9j The simulation produces a mobile 
glaciation with two distinctly different stages. Stage 1, the period of high average obliq- 
uity from 10 until 4 Ma ago, is characterised by ice thicknesses less than 400 m and rapid 
changes of the position where the maximum thickness occurs between all latitudes. By 
contrast, during stage 2, the period of low average obliquity from 4 Ma ago until today, the 
position of maximum thickness changes much less rapidly and flip-flops between the poles 
only (47% of the time at the north pole, 53% of the time at the south pole). The polar ice 
deposits grow almost monotonically to their present-day thicknesses, only interrupted by 
moderate shrinking around ~3.2, 1.9 and 0.7 Ma ago when maximum amplitudes of the 
main obliquity cycle of 125 ka occurred (see also Fig. [TJ. 

In order to illustrate this behaviour in more detail, Fig. [10] depicts the distribution 
of the simulated ice thickness for three selected time slices. The first time slice, 4.1 Ma 
ago (near the end of stage 1), is characterised by a low maximum ice thickness (74.4 m) 
which occurs in the low southern latitudes (at 6°S) and a continuous glaciation from the 
mid southern (47°S) to the low northern (22°N) latitudes. In addition, small south polar 
deposits extend from the pole to 83°S, whereas the north polar region is entirely ice-free. 

The second time slice, 0.61 Ma ago (in stage 2, towards the end of the most recent 
period of large obliquity amplitudes), shows polar ice deposits only moderately smaller 
than their present-day counterparts. However, the striking difference compared to the 
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Figure 10: (a) Simulation #6: Ice thickness H for three selected time slices (which cor- 
respond to the marks in Fig. [9^,.) The bottom panel shows thickness classes in order to 
highlight deposits of thin ice not visible in the upper panel. The class " <0 " means ground 
ice (see Sect. 2.7). (b) Observational data of the ice thickness of the present-day PLDs. 



They have been obtained by subtracting the MOLA surface topography ( Zuber et al.||1998 



Smith et al. 1999 ) from the basal topography computed by a smooth extrapolation of the 



ice-free ground (Greve et al. 2004) and subsequent zonal averaging. 



present is the existence of continuous, at least metres-thick ice deposits equatorward to 
38°N and 52°S, respectively. These mid-latitude deposits are very mobile; merely 20 ka 
earlier (0.63 Ma ago), they are entirely absent in the northern hemisphere while extending 
equatorward to 24° S in the southern hemisphere (not shown). 

The third time slice is the present (t = 0), with an obliquity of 25.2° following a 
~ 0.3 Ma period with only small changes (within less than 5°). As already discussed above, 
for the present, the simulation predicts bipolar ice deposits which match the observed 
volumes within ~ 3% and the ice thicknesses at the poles within ~ 20%. The bulk of the 
deposits with thicknesses > 100 m extend to 80°N and 77°S, which also agrees well with the 
observed values of 80°N and 73°S, respectively. Metres-thick deposits follow equatorward 
to 57°N and 70°S. This is still compatible with the reality in the south, while it contradicts 
the reality in the north where such deposits are not observed. 



4 Discussion and conclusion 



The simulations with constant orbital parameters presented in Sect. 3.1 confirm the in- 



tuitive idea that low obliquities favour deposition of water ice in high latitudes and vice 
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versa. An interesting additional finding is that the polar ice deposits for relatively low 
obliquities can either occur at one pole only or at both poles, depending on the asymmetry 
of the seasons in the two hemispheres. 

The more realistic best-fit simulation #6 of Sect. 32 with time-dependent orbital forc- 
ing from 10 Ma ago until today produces a very good agreement with the observed present- 
day PLDs. It predicts a mobile glaciation with two distinct stages. During stage 1, from 
10 to 4 Ma ago, ice thicknesses never extend 400 m, and ice is readily exchanged be- 
tween all latitudes. This exchange is mainly controlled by obliquity, polar deposits being 
again favoured by relatively low obliquities and lower-latitude deposits by peak obliquities. 
During stage 2, from 4 Ma ago until today, the north and south polar ice deposits grow 
essentially monotonically and reach their maximum thicknesses at the present. In partic- 
ular during the three periods of large obliquity amplitudes around ~3.2, 1.9 and 0.7 Ma 
ago, the polar deposits continue in thin, very mobile ice deposits which extend far into the 



mid latitudes at times. The latter result agrees with the findings by Head et al. (2003) 



who report evidence for "ice ages" on Mars during the period from about 2.1 to 0.4 Ma 
ago when the obliquity regularly exceeded 30°. According to the authors, the conditions 
during this period favoured the deposition of metres-thick, dusty, water-ice-rich material 
down to latitudes of ~ 30° in both hemispheres. 

A limitation of the results of this study must be noted. The estimated surface ages of 



the northern (at most 0.1 Ma) and southern PLDs (about 10 Ma) by Herkenhoff and Plaut 



(2000), which are based on crater statistics, are consistent with the simulated growth of 
ice deposits during the last 4 Ma for the north, but not for the south polar region (see 
Fig. 10). This may be related to the fact that most of the southern PLDs are covered by 
dust, whereas the water ice of their northern counterpart is exposed to the atmosphere. 
Consequently, at least for the present-day situation, the northern PLDs can readily ex- 
change water with the atmosphere, whereas the exchange is blocked to an unknown extent 
for the southern PLDs. This problem requires further attention. Nevertheless, we con- 
clude that the model MAIC-2 is a very useful tool which, despite its simplicity, can provide 
substantial insight in the evolution of the Martian surface glaciation over climatological 
time scales. 
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A Discrete formulation 

A.l Numerical grid 

MAIC-2 is a spatially one-dimensional model which features a dependence on latitude only. The 
grid points are located at monotonically increasing latitudes, 

<Pl, l = 0,...,L, (25) 
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where tpo = —ir/2 (south pole) and ifL = ir/2 (north pole). The generally non-equidistant spacing 
between subsequent grid points is 

Atpt = ipi - ipi-i , l = l,...,L. (26) 

Further, we define the latitudes in between grid points (at cell boundaries) as 

fl±i/2 = g ' ( ' 

Time is discretised uniformly by 

t n = t° + nAt, n = 0,...,N, (28) 
where t° is the initial time, At the time step and t = t° + N At the final time. 

A. 2 Surface temperature, evaporation, condensation 

Numerical evaluation of the LIT [Eqs. ([6])-(10)] and evaporation [Eqs. (12)-(17)] schemes is 
essentially straightforward and need not be detailed. The daily cycle of the surface temperature 



[Eq. (15)] for computing evaporation rates is sampled with the sub-daily time step At~oc = g SOi - 



As for condensation, the discrete version of the condition ( 18 ) is 

1 / n P sa ,t(Tp) 



= { At \ 1 g J (29) 

, otherwise . 



A. 3 Instantaneous mixing 



As mentioned in Sect. 2.6, instantaneous mixing of water vapour in the atmosphere can be 
described as the limit of infinite diffusivity, K — > oo. A finite- volume discretisation of the dif- 
fusion equation with finite diffusivity is described in an earlier version of this paper (archived 
at arXiv:0903. 2688^1 [physics.geo-ph]); however, the limit K — > oo cannot be carried out nu- 



merically with that scheme. Instead, a two-step procedure has been devised for the case of 
instantaneous mixing. 

In the first step, for any point of the model domain (I = 0, . . . , L), let us compute predictors 
of the water content at the new time step, uf, by ignoring the water transport, 

^ ~ f ~ 1 = - Cf . (30) 

In the second step, the resulting total amount of water shall be mixed over the entire planet. 
The total amount of water is obtained by integrating over the grid cells and summing up, 

¥>l/2 L l fl + 1/2 

2ttR 2 u}q / cos ipd(p + 22 27ri? 2 w[ l / cos ip dip 
J ... J 

PL 



^tot 



+ 2itR 2 uj1 J cos <p dip 



PL-l/2 

2ttR 2 < cb% (1 + sinv? 1/2 ) + ^ cbf (sin^ +1/2 - sin^_ 1/2 ) 
I i=i 

+ w2(l-sin^ L _ 1/2 )} , (31) 
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where R denotes the radius of the planet (R = 3396 km). For all points / = 0, . . . ,L, the new 
water content follows by division by the surface of the planet, 

w ^ ( L-l 



A. 4 Ice evolution 



i=i 



+ u n L {l-sm<p L _ l/2 )) . (32) 



The discretisation of the ice-thickness equation (22) is straightforward. By using an Euler back- 
ward step for the time derivative, we obtain 

Tjn Tjn—1 /~in rn 

At Pice 
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